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Abstract 

Red colobus monkeys, due to their sensitivity to environmental change, are 
indicator species of the overall health of their tropical rainforest habitats. As a 
result of habitat loss and overhunting, they are among the most endangered 
primates in the world, with very few viable populations remaining. Tradition- 
ally, extant indicator species have been used to signify the conditions of their 
current habitats, but they have also been employed to track past environmental 
conditions by detecting previous population fluctuations. Kibale National Park 
(KNP) in Uganda harbors the only remaining unthreatened large population of 
red colobus. We used microsatellite DNA to evaluate the historical demography 
of these red colobus and, therefore, the long-term stability of their habitat. We 
find that the red colobus population throughout KNP has been stable for at 
least -40,000 years. We interpret this result as evidence of long-term forest sta- 
bility because a change in the available habitat or population movement would 
have elicited a corresponding change in population size. We conclude that the 
forest of what is now Kibale National Park may have served as a Late Pleisto- 
cene refuge for many East African species. 



Introduction 

During the Pleistocene and into the Holocene, there have 
been large faunal and floral changes due both to anthro- 
pogenic factors and climate change (Gill et al. 2009; 
Nogues-Bravo et al. 2010; McGlone 2012; Prescott et al. 
2012; Rule et al. 2012). Understanding the magnitude of 
each of these issues has been incredibly challenging. One 
method to detect environmental change is to determine 
when species populations have fluctuated (Goossens et al. 
2006; Schneider et al. 2010; deBruyn et al. 2011; Ting 
et al. 2012). Although species may respond differently to 
environmental changes (Arenas et al. 2012; Lorenzen 
et al. 2011), indicator species are particularly sensitive to 
minor changes in environmental condition (Lindenmayer 
et al. 2000), and in many cases, it is understood how they 
will respond to environmental changes. In fact, these spe- 
cies have been important for monitoring the status of our 



natural environment and even the presence of other spe- 
cies. Although, traditionally, extant indicator species have 
been used to indicate the conditions of their current 
environment, by looking at past population demographics 
(e.g., population bottlenecks or expansions), they may 
also be used to indicate the history of their habitat. 

Red colobus monkeys, genus Procolobus[Piliocolobus] 
(Oates and Davies 1994; Grubb et al. 2003), are indicator 
species (Chapman and Lambert 2000; Waltert et al. 2002; 
Struhsaker 2005, 2010) because they are sensitive to habi- 
tat change (Oates and Davies 1994; Waltert et al. 2002). 
Their populations are thought to signify the overall con- 
dition of their tropical rainforest. With rapid habitat loss 
across equatorial Africa and hunting in many areas, they 
are considered one of the most endangered primates in 
the world (Struhsaker 2005; Mittermeier et al. 2009). 
Although the taxonomy of red colobus has been 
under considerable debate for some time, there are 18 
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morphotypes currently recognized across equatorial Africa 
(Grubb et al. 2003; Ting 2008). In 2008, 13 morphotypes 
were evaluated by the IUCN and all but one was in 
decline and many were critically endangered. Unfortu- 
nately, little is known of the natural history of many of 
these morphotypes and the majority of red colobus habi- 
tat remains unprotected making the future of this genus 
very questionable (Struhsaker 2005). 

Procolobus [Piliocolobus] rufomitratus is found in Kibale 
National Park (KNP) in Western Uganda (Fig. 1), and is 
one of the few large groups from the entire distribution 
of this genus that is presently not threatened. The esti- 
mated census size (N c ) of the red colobus in KNP is 
-17,000 individuals (Struhsaker 2005), which is indicative 
of a healthy tropical rainforest that occupies a relatively 
small area of -795 km 2 . Kibale National Park supports an 
extremely high level of biological diversity with more than 
351 species of trees, 335 species of birds, and an unprece- 
dented 13 species of primates (Uganda Wildlife Authority 
2011). Although there has been much deforestation in 
Uganda (Howard et al. 2000), KNP has remained rela- 
tively stable since early monitoring began in the 1970s 
(Struhsaker 1975; Southworth et al. 2010; Naughton-Treves 
et al. 2011). This upland tropical rainforest was desig- 
nated a national park in 1993 to preserve its biological 
diversity. The park is a major destination for eco-tourists 
and has been at the forefront of research on African forest 
primates for over 40 years (Chapman et al. 2005, 2010a,b; 
Box et al. 2008; Struhsaker 2010). However, we have little 
information about the history of Kibale forest, particu- 
larly during the glacial and interglacial periods of the Late 
Pleistocene and into the early Holocene, where we know 
surrounding areas were affected (Scholz et al. 2007; Stager 
and Johnson 2008). Similarly, little is known about 
anthropogenic affects of the forest, apart from the arrival 
of the Bantu-speaking people around -2,300 years ago 




Figure 1. Red colobus monkey, Procolobus [Piliocolobus] 
rufomitratus, from Kibale National Park, Uganda. 
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who are thought to have caused major changes in the for- 
est habitat surrounding Kibale forest (Taylor et al. 1999). 
As red colobus are very sensitive to habitat change, 
evidence of past population size changes would likely 
indicate changes in Kibale forest and shed light on past 
factors that affected this area. 

Microsatellites (DNA regions of short tandem repeats) 
provide a rich source of polymorphic nuclear DNA mark- 
ers for studying the genetic structure and diversity of a 
population's gene pool (Jarne and Lagoda 1996; Cavalli- 
Sforza 1998; Ellegren 2004). Microsatellite data can be 
used to look for evidence of past population decline or 
growth by applying methods based on the coalescent the- 
ory and can provide unique insights into the long-term 
history of a population. Currently, it is not known if all 
of the red colobus in KNP form one large population 
made up of many interbreeding groups or if there are dif- 
ferent populations within the park. This information 
gained from microsatellite data for the red colobus in 
KNP can be compared to other groups of red colobus 
that are currently declining. 

We used microsatellites from six groups of red colobus 
living throughout KNP to first estimate the number of 
populations in the park using three different analytical 
approaches, including Bayesian model testing. We then 
examined the KNP population for evidence of historical 
fluctuations in effective population size (N e ) using three 
complementary approaches that included Extended Bayes- 
ian Skyline Plots (EBSP) (Heled and Drummond 2008). 
These models ranged from the familiar stepwise mutation 
model (SMM) to more complex ones that account for 
various properties of microsatellite mutation (Sainudiin 
et al. 2004). Finally, we integrated our population genetic 
results with field and paleoenvironmental data to infer 
the history of the forest of Kibale National Park during 
the Late Pleistocene through the Holocene. 

Materials and Methods 
Sample collection 

Blood and/or fecal material was obtained from 64 indi- 
viduals of the Large Mikana and Small Camp groups, 
which are both part of a long-term field study that was 
started in 2003. Two to four fecal samples were collected 
from each individual that was identified by distinguishing 
characteristics (e.g., scars, broken tails, hair color, and 
birthmarks). In an effort to sample from throughout the 
park, fecal specimens were also collected from 21 individ- 
uals from the Dura (N = 5), K30 (N = 5), Mainaro 
(N = 5), and Sebatoli (N = 6) groups (Fig. 2). The blood 
samples were stored on classic FTA cards, and the fecal 
specimens were treated with the two-step procedure 
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Figure 2. (a) Distribution of red colobus across north central Africa (black shading; Ting 2008), with the red star highlighting Kibale National Park 
(KNP). (d) Uganda with the location of KNP in red. (c) KNP and the locations of the six groups on the days of their sampling along with their 
sample sizes (modified from a map by A. Jacob). 



described in Nsubuga et al. (2004). Total genomic DNA 
was extracted from the blood and fecal specimens with 
the FTA protocol and Qiagen QIAmp" Stool Kit (cat: 
51504) according to their manufacturers' recommenda- 
tions. To obtain adequate DNA for genotyping, DNA was 
extracted from two to four fecal specimens per individual. 

Eleven microsatellite loci from the human genome, 
found to be polymorphic in red colobus (A.B. Di Fiore, 
pers. comm.), were selected for genotyping (Table 1). 
These loci were amplified with a three-primer nested PCR 
approach (Schuelke 2000), run on an AB3730xl automated 
sequencer (96 capillary), and scored with GeneMarker® 
(Softgenetics, State College, PA). All loci were amplified 
multiple times from multiple extracts per individual to 
ensure accuracy of genotypes. An earlier version of the 
red colobus dataset was used by Wu and Drummond 
(2011) to test the performance of their EBSP models for 
microsatellites. This earlier dataset included only 31 of 



the 85 red colobus, which are now represented in the 
current matrix. Our present EBSP results supersede these 
previous ones as they are based on all 85 individuals 
and on direct comprehensive tests for the historical 
demography of red colobus in KNP (Appendix Al). 

Standard statistics 

Although the six groups were initially treated as separate 
for the calculation of summary statistics, they were ulti- 
mately combined according to our series of four tests for 
population number and ongoing field observations with 
red colobus at KNP (see below). The combined dataset 
was checked for genotyping errors with MICRO- 
CHECKER, version 2.2.3 (van Oosterhout et al. 2004). 
Estimates of within-group diversity [observed and effec- 
tive numbers of alleles (n a and n e ) and observed and 
expected heterozygosity {H 0 and H e ) ] were calculated with 
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Table 1. Standard statistics for the 10 microsatellite loci of the 85 red colobus (Procolobus [Piliocolobus] rufomitratus) from Kibale National Park 



in Uganda. 


Locus 3 




n a 


r/ P 


H Q 


He 


p i e 


D14S306 


84 


6 


4.25 


0.76 


0.77 


0.01 


D3S1766 


82 


9 


6.31 


0.76 


0.85 


0.11 


D2S1399 


84 


10 


6.57 


0.91 


0.85 


-0.06 


D7S1817 


78 


9 


4.88 


0.69 


0.8 


0.14 


D20S206 C 


79 


7 


5.08 


0.79 


0.81 


0.03 


D8S60 


82 


4 


1.91 


0.54 


0.48 


-0.12 


D8S165 


85 


3 


2.07 


0.55 


0.52 


-0.06 


D1S207 


85 


14 


6.12 


0.74 


0.84 


0.12 


C2A 


85 


7 


2.39 


0.55 


0.59 


0.06 


D561457 


82 


8 


3.62 


0.67 


0.73 


0.08 


Average d 


82.6 


7.70 


4.32 


0.70 


0.72 


0.036 




(78-85) 


(3-14) 


(1.91-6.57) 


(0.54-0.91) 


(0.48-0.85) 


(-0.017-0.082) 



a Locus nomenclature follows that of the human genome. 

b Standard statistics include: (W): number of genotyped individuals per locus, (n a ) observed number of alleles, (n e ) effective number of alleles, (H 0 ) 
observed heterozygosity, (H e ) expected heterozygosity, and (F /s ) Wright's F ts . 

This locus is the only one that deviates significantly from Hardy-Weinberg expectations after Bonferroni correction. 

d The average for F ts corresponds to its weighted mean and 95% bootstrap confidence interval across loci. The averages for the other statistics 
refer to their arithmetic means and ranges. 

e F B estimates (means and 95% confidence intervals, in parentheses) are 0.10 (0.02-0.17) for Large Mikana and -0.01 (-0.07-0.04) for Small 
Camp. Thus, total F, s is not significantly negative for either group with a sample size of >6 individuals. A significantly negative F /s is indicative of a 
heterozygote excess and may be related to female philopatry (Chesser 1991). However, dispersal in red colobus is female biased (Struhsaker 
2010) and so these results are not unexpected. Conversely, the significantly positive F s for Large Mikana is indicative of a heterogyzote deficit 
and may be related to locus-specific underdominant selection (Miyamoto et al., unpubl. manuscript). 



POPGENE, version 1.3.1 (Yeh and Boyle 1997). Each 
locus was tested against its Hardy-Weinberg expectations 
with the G-test of POPGENE. Each pair of loci was 
assessed for genotypic disequilibrium with the G-test and 
4,500 permutations in FSTAT, version 2.9.3 (Goudet 
2002). Wright's F IS was estimated for each locus with 
10,000 randomizations using FSTAT and across all loci 
with 10,000 bootstrap replications in GDA, version 1.1 
(Lewis and Zaykin 2002). To correct for multiple compar- 
isons, a sequential Bonferroni correction was applied in 
all sample-wide series of tests (Rice 1989). 

Population number 

Three tests were conducted to determine whether the six 
groups make up a single interbreeding population or 
multiple separate populations. The first approach relied 
on estimates of R ST and 4Nm as calculated with RST- 
CALC, version 2.2 (Goodman 1997) for all six groups 
and between all 15 pairs of groups. The significance of 
the R S T estimates was determined with 10,000 permuta- 
tions. The second approach used STRUCTURE, version 
2.3.3 (Pritchard et al. 2000) and an admixture model with 
correlated allele frequencies among populations. The 
number of populations (K) was varied from one to six 
(one population per group). Each a priori setting of K 
was run 10 independent times with each run consisting of 



1 x 10 s generations and a burn-in of 1 x 10 7 . Results 
were then summarized with STRUCTURE HARVESTER 
version 0.6.92 (Evanno et al. 2005). 

Finally, we used log Bayes Factor (In BF) testing with 
MIGRATE version 3.1.6 (Redelings and Suchard 2005; 
Beerli and Palczewski 2010). We used a sequential hierarchical 
approach for the testing of simpler to more complex 
models with In BF. If the In BF test for the current 
simpler model was decisive, we stopped there; otherwise, 
we would continue with a comparison of the more com- 
plex model to one with an extra population, O, or 4Nm. 
For this analysis, we began with a comparison of the sim- 
plest one-population model (with its one free parameter, 
O) to its closest competitor with two free parameters (the 
two-population model with mean O and average 4Nm). 
The two groups for the two-population model were 
chosen according to the longest internal branch of their 
neighbor-joining tree as generated with the A/j 2 pairwise 
distances of the six groups (Goldstein et al. 1995). The 
two-population model assumed a single mean Q for both 
groups and symmetrical migration between them (i.e., an 
average 4Nm). The calculation of log marginal likelihoods 
(In mL) for the two models was done using thermody- 
namic integration and its Bezier approximation (Beerli 
and Palczewski 2010). These calculations relied on static 
heating with the temperatures of the three heated chains 
set at 1.5, 3.0, and 100,000. Each hot and cold chain was 
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run for 1 x 10 s generations per locus with samples taken 
from the cold chain every 2,500 steps after a 10% burn- 
in. Uniform priors of 0-50 and 25-400 were used for 0 
and 4Nm, respectively. Three independent runs under 
both the one- and two-population models were per- 
formed to ensure the convergence of the final results. All 
MIGRATE analyses were done using the Brownian 
approximation of the SMM. Rate variation among the 
loci was accounted for by the incorporation of mutation 
rate modifiers that were calculated according to the num- 
ber of observed alleles. 

Population size changes 

We performed four tests to check for historical changes 
in population size. The first approach relied on BOT- 
TLENECK, version 1.2.02 using the one-tailed Wilcoxon 
test for heterozygote excess and the module shift descrip- 
tor of the allele frequency distribution (Luikart et al. 
1998). The one-tailed test was used for the former to 
increase the power of detecting a heterozygote excess, 
which is indicative of a recent population bottleneck. 
The underlying distribution of H 0 in the Wilcoxon test 
was generated with 10,000 simulations under the two- 
phase mutation (TPM) model of microsatellite evolution. 
In this model, most mutations change the allele length 
by a single step, whereas a small proportion changes the 
allele size by two or more steps (Di Rienzo et al. 1998). 
The variance among multiple steps in the BOTTLENECK 
simulations was set to 12 as recommended by Piry et al. 
(1999). Then, M_P_VAL (Garza and Williamson 2001) 
was also used to test for a recent bottleneck. We tested 
the significance of the M ratio with 10,000 simulations 
under TPM in M_P_VAL. The proportion of single-step 
mutations and mean step size in the M_P_VAL simula- 
tions were set to 0.88 and 2.8, respectively, following the 
empirical recommendations of the program's authors. In 
the BOTTLENECK and M_P_VAL analyses, "recent" 
refers to a bottleneck where (e.g.) the population has 
decreased within the last 250 generations (i.e., 
1,250 years for red colobus given a generation time of 
5 years; see below) to N e of 50 (Cornuet and Luikart 
1996; Luikart et al. 1998). The recommended settings of 
Piry et al. (1999) and Garza and Williamson (2001) were 
adopted for the BOTTLENECK and M_P_VAL simula- 
tions, because they are empirically based on a number of 
different studies. 

The second approach used MSVAR, version 1.3 to 
estimate the starting time (f) for a change in the ancestral 
(Nj) to current (N 0 ) population sizes (Beaumont 1999). 
The ratio of N 0 to N, (Ng/Nj or r) was also calculated 
from their Markov chain Monte Carlo (MCMC) samples. 
The MSVAR analyses were performed with the exponen- 
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tial model of population growth and a generation time of 
5 years for red colobus as estimated from their ages at 
sexual maturity (Struhsaker pers com.). Each analysis was 
run for 2 x 10 10 generations with samples taken every 
500,000 steps after a 10% burn-in. The log 10 means and 
variances for the hyperprior means of N 0 , N lt and t were 
set to (3.6, 1.0), (3.6, 1.0), and (5.0, 1.0), respectively 
(other prior and hyperprior settings were also tried; see 
below). The log 10 means and variances for their hyperpri- 
or variances were all set to (0.0, 0.5). Three independent 
runs were performed to ensure convergence. 

The final approach estimated N e through time using 
EBSP (Heled and Drummond 2008) in BEAST version 
1.5.4. The Extended Bayesian Skyline Plot is a non-para- 
metric piecewise model, which estimates the population 
function directly from the data instead of fitting a para- 
metric model (e.g., exponential growth) specified a priori. 
Unlike the original Bayesian Skyline Plot (Drummond 
et al. 2005), EBSP combines information from multiple 
loci and estimates the number of population size change 
points, which is specified by the user as a fixed value in 
the Bayesian Skyline Plot. The analysis was performed 
with 12 microsatellite models (Sainudiin et al. 2004; 
Table 2), with four independent runs carried out for each 
implemented model (Wu and Drummond 2011). The 
observed maximum length of any allele was 33 repeats, so 
we assumed an upper bound of 35. The average mutation 
rate of all loci was fixed to one, but the relative rates were 
estimated to account for rate variation across the loci. At 
time point tp 9j was assumed to come from an exponen- 
tial distribution with an unknown mean parameter (4)), 
to which a one-on-x prior was applied. A Poisson prior 
with mean X = 0.6931 was used for the number of 
changes (t/0 in the population history (Heled and Drum- 
mond 2008). Each analysis was based on 8 x 10 s generations 
with samples taken every 4 x 10 steps after an appropri- 
ate burn-in of 10-15%. 

We next completed BEAST analyses of the microsatellite 
data, assuming a model of constant population size. These 
constant population size analyses using BEAST were com- 
pared to the EBSP runs to test further the hypothesis of a 
stable population of red colobus at KNP. These additional 
BEAST analyses were conducted with the best-fit microsat- 
ellite model from the EBSP comparisons (PU2; see Table 2) 
as well as with the EU1 model, which best approximates 
the SMM. We applied coalescent priors to all sampled 
genealogies of the microsatellites, which assume a constant 
population size that is shared across all loci. The common 
population size parameter was given the one-on-x prior, 
f_X(x) = 1/x, which is the Jeffrey's prior for a constant 
population size. Log BF were then used to compare the fits 
of the constant population size model to those of EBSP as 
described in the next paragraph. 
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Table 2. Twelve models of microsatellite evolution. 



Mutation 
Model 3 


Ln marginal likelihood 
(Ln mL) b 


Ln Bayes factor 0 (Ln BF) 


Mean coalescent time (7" COAL ) d 


Model description 3 


PU2 e 


-942.14 


0 


AO AAA 1 '~>~l AAA 

98,000 (3/, 000- 


1 — T f AAA\ 

-1 /6,000) 


Proportional-rate, unbiased, two-step 


EC2 


-946.28 


4.15 


105,000 (40,000- 


-183,000) 


Equal-rate, constant-bias, two-step 


EL2 


-950.43 


8.29 


121,000 (51,000- 


-212,000) 


Equal-rate, linear-bias, two-step 


PC 2 


-951.90 


9.77 


114,000 (41,000- 


-214,000) 


Proportional-rate, constant-bias, two-step 


PU1 


— 957 99 


1 5 86 


153,000 (61,000- 


-267,000) 


Pmnrii'+irinal-i'3+o i mhi^coH nno-cton 
r I (J|jui uui i a I I a Ltr, u I luiaitru, ui itr 2>itr|J 


EU2 


-969.63 


27.50 


107,000 (42,000- 


-182,000) 


Equal-rate, unbiased, two-step 


PC 1 


-980.49 


38.36 


163,000 (62,000- 


-286,000) 


Proportional-rate, constant-bias, one-step 


EU1 


-987.95 


45.81 


1 58,000 (69,000- 


-277,000) 


Equal-rate, unbiased, one-step 


EC1 


-992.24 


50.11 


154,000 (70,000- 


-257,000) 


Equal-rate, constant-bias, one-step 


PL2 


-997.90 


55.77 


129,000 (51,000- 


-235,000) 


Proportional-rate, linear-bias, two-step 


PL1 


-1006.89 


64.76 


197,000 (71,000- 


-349,000) 


Proportional-rate, linear-bias, one-step 


EL1 


-1023.05 


80.92 


180,000 (76,000- 


-310,000) 


Equal-rate, linear-bias, one-step 



a AII models have been described in Sainudiin et al. (2004), with the exception that we use a simple logistic regression (rather than simple linear 
regression) to model a directional bias of mutation (C and L below). These 12 models vary (1) according to whether the mutation rate is indepen- 
dent of allele length (E = equal, rate is independent of length and P = proportional, rate is proportional to microsatellite length); (2) according to 
whether the probability of a contraction is equal to that of an expansion (U = unbiased, equal probabilities for the two, C = constant bias, the 
probabilities of contraction and expansion are not equal, but remain constant and are independent of allele length, and L = linear bias, unequal 
probabilities for the two events, which now depend on allele length), and (3) according to whether a mutation can change the length of an allele 
by more than one repeat (1 = single step, only mutations of a single repeat are allowed and 2 = multi-step, a mutation can change the length of 
an allele by one or more repeats). Thus, the PU2 model accounts for a microsatellite mutation rate that varies according to allele length, equal 
rates of contraction and expansion, and mutations of multiple, as well as single, repeats. 
b The 12 models are listed from best to worst according to their In mL. 
Relative to the best (PU2) model. 
d 95% HPD in parentheses. 

The PU2 model receives >98% of the relative probability (i.e., the total marginal likelihood for all 12 models). 



Log Bayes Factors between each pair of microsatellite 
mutation and/or population (constant population size vs. 
EBSP) models were computed by TRACER version 1.4.1 
(Rambaut and Drummond 2007) for model selection. 
Here, the harmonic means of the sampled log likelihoods 
(ln L) were used to calculate the ln mL of the models. 
The models were then ranked from best to worse accord- 
ing to their ln BF. 

The MCMC sampling in each run with STRUCTURE, 
MIGRATE, MSVAR, and BEAST was monitored in TRA- 
CER to confirm its mixing and convergence. Different run 
conditions, priors, and/or hyperpriors were also examined 
in addition to those reported above for the final produc- 
tion runs. These experiments were performed to establish 
the final settings of the production runs (e.g., lengths, 
burn-in, and underlying distributions) and to verify that 
small changes did not result in significantly different esti- 
mates. The log files for the final production runs and for 
these additional experiments are available upon request 
from the corresponding author. Readers can inspect these 
log files with TRACER to further assess the posterior 
probabilities, mixing, and convergence of the different 
population, substitution, and MCMC parameters. 

For these analyses, we used the mutation rate, 
5 x 10~ 4 , calculated from human data using detailed 



microsatellite mutation models (Whittaker et al. 2003). 
However, another mutation rate, 1 x 10~ 4 , has also been 
used in microsatellite studies of the population genetics 
of humans and other primates (Ruiz-Garcia 2005; Payseur 
et al. 2011). Changing our mutation rate would not 
change the overall findings of the study, but would 
increase the coalescent times (T C oal defined by the height 
of the tallest tree at each step of the MCMC) and N e by a 
factor of five. We chose to present our results using the 
Whittaker et al. (2003) rate because it supports smaller, 
and thus, more conservative estimates of central impor- 
tance to our conclusions and to conservation biology. For 
example, these results offer smaller estimates for the age 
of the red colobus population, and thereby Kibale forest, 
and for minimum viable population size. 

Results and discussion 
Basic statistics 

The original microsatellite dataset consisted of 11 loci 
from 85 red colobus living in KNP. Locus D17S1290 was 
the only locus with a significant F IS after Bonferroni cor- 
rection and it was significantly positive (indicative of a 
heterozygote deficit). This locus was also the only locus 
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with both a significant number of null alleles and a 
homozygote excess due to scoring errors according to 
MICRO-CHECKER. Due to these discrepancies, this locus 
was excluded from all further analyses, including the fol- 
lowing summary statistics (Table 1). An average of 7.70 
alleles was scored from each locus, and heterozygosities 
averaged 0.70 for H 0 and 0.72 for H e . None of the 45 
pairs of loci were in genotypic linkage disequilibrium 
after Bonferroni correction. In turn, only one locus 
{D20S206) deviated significantly from Hardy-Weinberg 
expectations, but this departure was not due to either a 
heterozygote deficit or excess according to its non-signifi- 
cant F IS . This result of a significant Hardy-Weinberg 
test, but insignificant F IS was due to the existence of a 
heterozygote excess for the genotypes of some alleles, but 
a heterozygote deficit for others. 

Number of populations 

All three approaches for examining the number of popu- 
lations of red colobus from KNP supported a single 
population. First, the Rst and 4Nm estimates across all 
six groups were negative according to RSTCALC (—0.01 
and —25.28, respectively), with the former not deviating 
significantly from zero according to its permutation test 
(P > 0.10). None of the _R sr for the 15 pairwise compari- 
sons between groups were significant (after Bonferroni 
correction). Thus, the RSTCALC results reveal no 
significant genetic differences among the six groups, and 
indicate that the six groups can be considered as members 
of one panmictic population. 

Second, mean In L followed a U-shaped pattern as K 
(number of populations) was increased from 1 to 6 in 
the STRUCTURE analysis (Fig. 3). This pattern of ele- 
vated mean In L for larger values of K has been reported 
by the authors of STRUCTURE (Pritchard and Wen 
2003) as well as by others (Evanno et al. 2005). It has 
been attributed to inflated support for values of K greater 
than the true number of populations. Thus, of greater 
importance is the failure of the STRUCTURE analysis 
with K = 6 to clearly assign the 85 individuals of red col- 
obus to separate populations (Fig. 3). This failure, cou- 
pled with the earlier warnings of others, indicates that 
one population is more likely than two or more for the 
red colobus of KNP. 

Finally, the longest internal branch in the unrooted 
neighbor-joining tree with Afi 2 distances divided the six 
groups into two larger groups of Dura, K30, and Large 
Mikana versus Mainaro, Sebatoli, and Small Camp (see 
Fig. 2 for group names). The In mL for the one-popula- 
tion model was -8,196.98 according to MIGRATE. 
Conversely, the In mL for the two-population model (with 
its two groups as defined above by the neighbor-joining 
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Figure 3. (Above) Mean In L (circles) and their standard deviations 
(vertical lines) for K = 1 to 6 populations as obtained with 10 
independent runs per K with STRUCTURE. (Below) Relative 
probabilities of assigning the 85 individuals of red colobus to six 
(differently colored) populations. These STRUCTURE results clearly 
show that little to no support exists for the assignment of the 85 
individuals to separate distinct populations. 

tree) was -26,261.91. Thus, the two models differed by a 
In BF of 18,064.93. In the statistics literature, a In BF of >5 
is generally considered to be "decisive" for one model over 
another (Kass and Raftery 1995). Following this literature, 
the one-population model is thereby taken as strongly 
favored over the two-population model with mean O and 
symmetrical migration. Given this selection of the simpler 
one-population model, we stopped there. Otherwise, we 
would have continued with an evaluation of a two-popu- 
lation/three parameter model (e.g., one with asymmetrical 
migration and mean O). Taken together, the RSTCALC, 
STRUCTURE, and MIGRATE results suggest high levels 
of gene flow throughout the park indicating that the red 
colobus of KNP is one large interbreeding population. 
These results fit with observational data from the field. 
For example, we have observed three female dispersals and 
one male dispersal in a year between the Large Mikana 
and Small Camp groups. On the basis of the genetic and 
field evidence, all groups were analyzed as a single popula- 
tion in the calculation of the summary statistics (Table 1) 
and in the following analyses of population fluctuation. 
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Population fluctuations 

The results from BOTTLENECK and M_P_VAL showed 
little evidence of a recent bottleneck for the red colobus 
in KNP. The one-tailed Wilcoxon test for heterozygote 
excess was not significant (P = 0.410) and the allele fre- 
quency distribution was L-shaped according to BOTTLE- 
NECK. The M ratio was calculated as 0.839, which was 
not significant (P = 0.303) according to its simulations in 
M_P_VAL. This large M ratio agrees with the empirical 
estimates of other species with stable population sizes 
(Garza and Williamson 2001). 

The MSVAR results were less clear. The averages and 
95% highest posterior densities (HPD) for log 10 N 0 (cur- 
rent N e ), Ni (ancestral N e ), r (Nq/Nj), and t (starting time 
of the population size change) were (3.22, 2.60-3.75), 
(3.92, 3.24-4.67), (-0.70, -1.58-0.16), and (4.18, 2.54- 
5.97) according to the 36,000 MCMC samples of MSVAR. 
Thus, the 95% HPD for log 10 r included zero, which is 
the critical point for N 0 = Nj. This critical point fell 
within the right tail of the posterior distribution at P 
(two-tailed) = 0.092 (Fig. 4). Conversely, the transformed 
linear means and 95% HPD for r and t were (0.20, 0.03- 
1.45) and (15,136, 347-993,254 years ago), which suggests 
a decline in N e of -80% over the last -15,000 years. We 
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Figure 4. Distribution of log 10 r (log 10 Ng/N,) summarizing the 36,000 
MCMC samples from MSVAR. A/j represents the past population size 
and N 0 represents the size of the current population. The upper 
boundary of our 95% credibility interval overlaps log q0 r = 0 as shown 
by the vertical gray line. This is where W 0 = W,, which is indicative of no 
population size change over time. The log 10 r estimates of this 
histogram support a model of population decline over growth by a In BF 
of 20.518 (Storz & Beaumont 2002). Conversely, our study concludes 
in favor of a constant population given that Iog10 r = 0 falls within the 
95% credibility interval of these estimates. 



also note that Storz and Beaumont (2002) described a In 
BF test, which compares how well a model for population 
contraction fits the data relative to one for expansion. 
Nevertheless, because our 95% HPD for the log 10 
r included zero, we cannot reject the hypothesis of a con- 
stant population. Furthermore, the estimated mean and 
median for the start of this decline (-15,000 years ago) 
falls within a period when tropical forest regions were 
expanding as earth entered its current interglacial (Casta- 
neda et al. 2009). Thus, when coupled with the non-sig- 
nificant deviation of log 10 r from zero, we cannot exclude 
a constant population hypothesis. 

The EBSP for all models show a flat line for N e since 
mean T COAL (Fig. 5). The posterior probabilities of the 
demographic.popSizeChange was the highest when it was 
0 and all 95% credible intervals included 0 for all models 
suggesting no population size change in the past. Of 
these 12 models for microsatellite evolution, the propor- 
tional-length, two-step (PU2) model had the best fit to 
our dataset (Table 2). This model was supported over its 
closest competitor, the equal-length, constant-size, two- 
step (EC2) model, by a non-decisive In BF of 4.15 (Kass 
and Raftery 1995). The poorest model was the equal- 
rate, linear-size, single-step model (ELI; In BF = 80.92 
vs. PU2; Table 2). The mean T COAL was estimated 
between 98,000 and 197,000 years ago among the 12 
models, with the best model (PU2) predicting 
98,000 years ago. Similarly, MIGRATE, using the Brown- 
ian approximation of the SMM and the one-population 
model, calculated a mean Pcoal of 112,848 years ago 
(standard deviation = 71,414). The In mL for the con- 
stant population size analyses with the PU2 and EU1 
models were —943.44 and —988.68, respectively, and dif- 
fered from those of their EBSP counterparts with PU2 
and EU1 by In BF of 3.68 and 2.08, respectively 
(Table 2). These small In BF are non-decisive according 
to Kass and Raftery (1995). Thus, the near-equal fits of 
both constant population size models confirm that the 
EBSP graphs are indicative of a stable population for the 
red colobus in KNP. 

The EBSP approach is powerful and flexible, because it is 
not constrained a priori to a specific model of historical 
population size change (Heled and Drummond 2008). 
Thus, the EBSP recovery of a constant population for the 
red colobus in KNP is derived from the data rather than 
from the a priori selection of a particular growth model. 
This support for a constant population is based on the 
maximum clade credibility (MCC) genealogies of each mi- 
crosatellite locus (Appendix A2), which show a standard 
coalescent pattern of increasing waiting times between 
coalescences as one works backward in time (Wakeley 
2008). Thus, relatively few internodes of each MCC gene 
genealogy occur more than > 10,000 years ago. It is impor- 
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Figure 6. Significant positive correlation between H a and average 
allele length for seven tetranucleotide loci (Spearman Rank 
Correlation, P = 0.033). The average allele lengths for each locus are 
weighted by their allele frequencies. This comparison does not include 
the three dinucleotide loci of the full dataset. 
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Figure 5. Extended Bayesian Skyline Plots for the 10 microsatellite 
loci using the PU2 (proportional-length, unbiased, two-step) model (a) 
and EU1 (equal-rate, unbiased, one-step) model (£>). These skyline 
plots are extended to their mean coalescent times (T C oal) on the right 
end of the graph. The gray shading corresponds to the 95% HPD 
around the mean N B . PU2 exhibits the best fit to the data, whereas 
EU1 approximates the SMM (Table 2). Both models (as do the other 
ten) show similar skyline plots with flat lines that correspond to a 
stable population size over time coalescing between -100,000 and 
200,000 years ago. 



tant to recognize that this paucity of old internodes does 
not reflect a lack of support for a constant population with 
mean and median Tqoal of -100,000 years ago. Instead, 
this paucity of old internodes contributes to the signal for a 
constant population, rather than for some other model 
such as for an expanding population, where the coales- 
cences are expected to cluster closer to the base of each gene 
genealogy (Page and Holmes 1998; Hein et al. 2005). How- 
ever, in keeping with our conservative approach, we limit 



our estimate of the population stability in red colobus to at 
least ~40,000 years ago (i.e., to the lower boundary of the 
95% HPD of T C oal for the PU2 microsatellite model; 
Table 2). This conservative treatment of the duration of the 
population stability is also warranted given that T COAL is an 
estimate of the coalescent time of gene genealogies and not 
necessarily the most recent common ancestor of the popu- 
lation (Hein et al. 2005). 

The evolution of microsatellites is complex because it is 
characterized by various mutational biases (Goldstein and 
Pollock 1997). A common property is the positive rela- 
tionship between allele length and mutation rate, which is 
evident in our data (Fig. 6). Such complexities call for 
the development and use of more realistic evolutionary 
models for microsatellites, such as the PU2 model 
(Table 2). Nevertheless, the fact that even simple models 
(e.g., EU1 which approximates the SMM) produce skyline 
plots similar to those of complex models (Fig. 5) suggests 
that simpler models can also be of value, particularly 
when dealing with large microsatellite datasets that are 
computationally intensive. The maximum clade credibility 
trees scaled to years for the 10 loci for models PU2 and 
EU1 are in Appendix A2. 
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Effective and census population sizes 

Mean N e varied among the EBSP graphs for the 12 
micro satellite models from -2,500 for the two best models 
(PU2 and EC2) to -4,000 (Fig. 5). The BEAST analyses 
with a constant population size model supported mean 
N e of 2,517 (95% HPD of 1,347-3,710) and 3,677 (2,561- 
4,857) for the PU2 and EU1 models, respectively. In turn, 
the one-population analysis with MIGRATE (which also 
assumes a constant population size) supported mean N e 
of 4,055 (2,978-5,426). Collectively, these different analy- 
ses of population size agree on a mean N e of -2,500 to 
4,000 for the red colobus in Kibale. 

Using this range of 2,500-4,000, we estimate the NJ 
N c ratio to be -0.15 to 0.25 for this population. These 
numbers are lower than the previous NJN C estimate of 
-0.35 for the red colobus population in KNP (Struhsak- 
er and Pope 1991). However, this previous NJN C esti- 
mate is based on older census techniques and 
calculations of male and female reproductive success, 
and is currently thought to be an overestimate (T. 
Struhsaker pers com.). Importantly, a number of factors 
are thought to reduce N e with respect to census esti- 
mates: skewed sex ratio, overlapping generations, and 
variation in adult reproductive success (Frankham 1995). 
The red colobus at KNP are characterized by all of these 
demographic factors. They exhibit a skewed sex ratio of 
one breeding male to every 1.3 breeding females (Struh- 
saker and Pope 1991). Red colobus also have overlap- 
ping generations, their overall lifespan is thought to be 
around 9.3 years (Pope 1996) with an overall reproduc- 
tive lifespan of 8.8 years, and adults average one birth 
every 2 years (Struhsaker and Pope 1991). Combined, 
these factors likely explain the low NJN C ratio found in 
this population of red colobus. 

In conservation biology, a critical N e of -500 to 1,000 
is often cited as necessary for the genetic security and 
long-term survival of a population (Franklin and Frank- 
ham 1998). However, other studies have argued from the 
perspective of the nearly neutral theory that a minimal 
viable population size of -1,000 to 5,000 is more appro- 
priate (Lynch and Lande 1998). Our N e estimates for the 
red colobus in Kibale (-2,500 to 4,000) overlap with the 
latter values, and thereby suggest that these larger N e are 
important for the genetic integrity and long-term viability 
of populations, similar to other recent studies (Flather 
et al. 2011). 

History of red colobus and Kibale forest 

Red colobus (Procolobus[Piliocolobus]) are found through- 
out most of equatorial Africa in small, mostly declining 
populations (Fig. 2). The morphotype, P. [Piliocolobus] 
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rufomitratus belongs to a central/eastern African clade of 
red colobus (Struhsaker 1981b; Grubb et al. 2003) that 
radiated around 1.4 million years ago (Ting 2008). 
Procolobus [Piliocolobus] rufomitratus likely split from 
other Piliocolobus sometime after 600,000 years ago (Ting 
2008), and our results suggest that the KNP population 
of red colobus has remained stable for -40,000 years or 
more. 

Red colobus are an important indicator species of their 
habitat, such that evidence of a population size change 
can be interpreted as a change in their habitat (Struhsaker 
2005). Our results suggest that these monkeys have been 
in a forest approximately the same size as KNP during at 
least the last -40,000 years as changes in the available for- 
est environment, or population movement, would have 
led to a corresponding change in population size. Fur- 
thermore, in recent years, no population size change in 
red colobus has been detected from 1970 to 2006 suggest- 
ing the forest is at carrying capacity for the red colobus 
(Chapman et al. 2010a,b). Thus, the most parsimonious 
hypothesis is that the forest of Kibale National Park has 
been persistently occupied by this population of red colo- 
bus for the duration of their population stability. This 
suggests that the forest has served as a refuge for tropical 
rainforest species during the last -40,000 years and possi- 
bly even longer. 

The alternative hypothesis, that this population of P. 
rufomitratus occupied a forest other than KNP during 
this time and then moved en masse to this area, is less 
likely for several reasons. The lack of a dramatic drop in 
population size does not support a recent founder event 
or the possibility that these red colobus are the remnants 
of a larger and more widely dispersed population. If 
these events occurred, they must have done so earlier 
than -40,000 years ago. It is also unlikely that this popu- 
lation of red colobus migrated en masse from another 
forest similar in size and environment sometime during 
the Late Pleistocene without any significant change in 
effective population size (e.g., bottleneck or founder 
event). Therefore, we conclude that the population stabil- 
ity of red colobus in KNP during the Late Pleistocene 
reflects similar long-term stability of the Kibale forest 
environment. 

East Africa has experienced many paleoenvironmental 
and anthropogenic changes over the last 100,000 years. 
This geographic region has undergone many long periods 
of extreme drought during the Late Pleistocene, first from 
-135,000 to 70,000 years ago (Scholz et al. 2007), and 
again during the last glacial maximum -18,000 years ago 
when Lake Victoria completely dried up (Stager and 
Johnson 2008) and regional vegetation changed to scrub 
and montane grasslands (Jolly et al. 1997). During these 
periods, a major Central Refuge for tropical forest species 
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existed to the west of Kibale forest over the Rwenzori 
Mountains in the neighboring Democratic Republic of 
the Congo (Hamilton 1976; Jolly et al. 1997; Wronski 
and Hausdorf 2008). Similar periods of environmental 
variability existed during the Holocene and are evidenced 
by a population crash of Ugandan and Kenyan buffalo 
indicative of a drought -4,500 years ago (Heller et al. 
2008). There is further evidence of anthropogenic impacts 
to the forests surrounding Kibale forest starting with the 
arrival of the Bantu-speaking people -2,300 years ago and 
through subsequent shifts in their settlement patterns 
(Taylor et al. 1999). 

Despite this variability, there is evidence that upper 
montane forests may have persisted during these climatic 
cycles (Jolly et al. 1997). Critically, Kibale forest is considered 
a separate habitat patch that persists due to its favorable 
mid-elevation topography and humidity (Struhsaker 
1981a; Jolly et al. 1997; Moore 1998; Wronski and 
Hausdorf 2008). Our results for a constant red colobus 
population in KNP corroborate the hypothesis of a stable 
Kibale forest, despite the known environmental and 
human changes in East Africa during the Late Pleistocene 
and Holocene. 

Our results further suggest that a tropical rainforest 
patch of only -795 km 2 (i.e., the size of KNP) is sufficient 
to support a long-term viable population of red colobus. 
This conclusion is important given the rapid rate at 
which their habitat is disappearing and the lack of knowl- 
edge of what makes a suitably sized conservation area. 
However, this estimate of a minimal habitat size refers 
specifically to red colobus because other species differ in 
their ecology and behavior and will have other habitat 
requirements (e.g., it is known that a territorial top pred- 
ator such as a leopard needs a much larger habitat for the 
long-term survival of its population). 

We now call for critical tests of our hypotheses for 
an old and stable KNP population of red colobus and, 
thus, Kibale forest. Palynological studies from through- 
out KNP are needed to critically test the hypothesis of 
a continuous old tropical rainforest of approximately 
the same size as that currently found in the park (Jolly 
et al. 1997). There has been no known connectivity 
between the red colobus of KNP and any other group 
since the 1950s. However unlikely, population genetic 
studies are now needed to test whether gene flow is 
occurring or has occurred with these other populations 
in East Africa. These studies are particularly important 
given that tropical rainforest in East Africa was much 
more widespread in the past (e.g., during the Holocene 
Climatic Optimum -6,000 to 10,000 years ago; Jolly 
et al. 1997). They are also significant in light of recent 
simulation studies (e.g., Chikhi et al. 2010; Peter et al. 
2010), which show that immigration can result in a 
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false signal of a population decline (i.e., the -80% 
decrease at -15,000 years ago as suggested by MSVAR). 
Additional simulations are also now needed to validate 
the efficiency and robustness of EBSP to recover older 
population size changes with microsatellite data (Wu 
and Drummond 2011). Such critical tests and experi- 
ments are necessary for the final acceptance of our con- 
clusions for an old and stable population of red 
colobus and Kibale forest. 

Biodiversity is lost at an astonishing rate each year, 
particularly in the tropics where deforestation is rapidly 
extinguishing species habitats. In Uganda, it is estimated 
that 85% of the forests have disappeared (Howard et al. 
2000). Wise conservation efforts over the last 40 years in 
Kibale National Park have ensured that this area has per- 
sisted as a protected National Park (Box et al. 2008). Our 
research suggests that this forest has been stable for at 
least the last -40,000 years, during which it likely served 
as an important sanctuary for tropical rainforest species 
in East Africa. Furthermore, our findings indicate that 
species with effective population sizes of -2,500 to -4,000 
may be able to persist long term in stable environments, 
thereby providing hope for other species in decline or 
with patchy distributions. 

Kibale forest offers a rich source of information about 
Pleistocene refugia, which increases its importance as a 
major center for conservation and scientific research. 
With continued wise management, KNP can likely con- 
tinue to serve as a critical refuge for these species well 
into the future. 
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Supporting Information 

Additional Supporting Information may be found in the 
online version of this article: 

Appendix Al. Genotypes for 85 red colobus from six 
groups in Kibale National Park (KNP). Missing data are 
indicated by an X. Group locations in KNP are given in 
Figure 1; SM - Small Camp, LM - Large Mikana, K30 - 
K-30, DUR - Dura, SEB - Sebatoli, and MAIN - Mainaro. 
Appendix A2. The maximum clade credibility trees for 
each of the 10 loci for two models of molecular evolution, 
(I) our best-fit model - PU2 and (II) EU1 the model that 
most approximates a stepwise mutation model. The time 
scale is in terms of years. 

Please note: Wiley- Blackwell are not responsible for the 
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